Utilisation de plotdf

load("plotdf");
plotdf([b1*x-k1*x*y,k2*x*y-d2*y],[x,y],
[parameters,"b1=1,k1=1,d2=1,k2=1"],
[tstep,0.1],[nsteps,100],
[direction,forward],[x,0,10],[y,0,10])$

Euler implicite

UnPasEulerImplicite(x0,h):=block([x:x0[1],y:x0[2]],[
((h+1)*sqrt(h^2*y^2+(2*h^2*x-2*h^3+2*h)*y+h^2*x^2+(2*h^3-2*h)*x+h^4-2*h^2+1)
 +(-h^2-h)*y+(h^2+h)*x+h^3+h^2-h-1)
 /(h*sqrt(h^2*y^2+(2*h^2*x-2*h^3+2*h)*y+h^2*x^2+(2*h^3-2*h)*x+h^4-2*h^2+1)
  +h^2*y+h^2*x+h^3-h),
(sqrt(h^2*y^2+(2*h^2*x-2*h^3+2*h)*y+h^2*x^2+(2*h^3-2*h)*x+h^4-2*h^2+1)
 +h*y+h*x+h^2-1)
 /(2*h^2+2*h)]);
UPEI:UnPasEulerImplicite;

Cranck-Nicholson

UnPasCrankNicholson(x0,h):=block([x:x0[1],y:x0[2]],[
-((h*x-2*h-4)*sqrt(4*h^2*y^2+(8*h^2*x-4*h^3+16*h)*y+4*h^2*x^2+(4*h^3-16*h)*x
                            +h^4-8*h^2+16)
 +(2*h^2*x+4*h^2+8*h)*y+2*h^2*x^2+(h^3-4*h^2-12*h)*x-2*h^3-4*h^2+8*h+16)
 /(h*sqrt(4*h^2*y^2+(8*h^2*x-4*h^3+16*h)*y+4*h^2*x^2+(4*h^3-16*h)*x+h^4-8*h^2
                   +16)
  +2*h^2*y+2*h^2*x+h^3-4*h),
(sqrt(4*h^2*y^2+(8*h^2*x-4*h^3+16*h)*y+4*h^2*x^2+(4*h^3-16*h)*x+h^4-8*h^2+16)
 -h^2*y+2*h*x+h^2-4)
 /(h^2+2*h)
]);
UPCN:UnPasCrankNicholson;

Dormand et Prince 5-4

dopri54(x0,fct,h):=
block([a21:1/5,
       a31:3/40,a32:9/40,
       a41:44/45,a42:-56/15,a43:32/9,
       a51:19372/6561,a52:-25360/2187,a53:64448/6561,a54:-212/729,
       a61:9017/3168,a62:-355/33,a63:46732/5247,a64:49/176,a65:-5103/18656,
       a71:35/384,a72:0,a73:500/1113,a74:125/192,a75:-2187/6784,a76:11/84,
       b1:35/384,b2:0,b3:500/1113,b4:125/192,b5:-2187/6784,b6:11/84,b7:0,
       bc1:5179/57600,bc2:0,bc3:7571/16695,bc4:393/640,
           bc5:-92097/339200,bc6:187/2100,bc7:1/40,
       x01,x02,x03,x04,x05,x06,x07,
       f01,f02,f03,f04,f05,f06,f07,
       x1,xt1
       ],
x01:x0,
f01:fct(x01),
x02:x0+h*a21*f01,
f02:fct(x02),
x03:x0+h*(a31*f01+a32*f02),
f03:fct(x03),
x04:x0+h*(a41*f01+a42*f02+a43*f03),
f04:fct(x04),
x05:x0+h*(a51*f01+a52*f02+a53*f03+a54*f04),
f05:fct(x05),
x06:x0+h*(a61*f01+a62*f02+a63*f03+a64*f04+a65*f05),
f06:fct(x06),
x07:x0+h*(a71*f01+a72*f02+a73*f03+a74*f04+a75*f05+a76*f06),
f07:fct(x07),
x1:x0+h*(b1*f01+b2*f02+b3*f03+b4*f04+b5*f05+b6*f06+b7*f07),
xt1:x0+h*(bc1*f01+bc2*f02+bc3*f03+bc4*f04+bc5*f05+bc6*f06+bc7*f07),
[x1,xt1]
)$